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Abstract 

A new algorithm is introduced to integrate the equations of rotational 
motion. The algorithm is derived within a leapfrog framework and the 
quantities involved into the integration are mid-step angular momenta 
and on-step orientational positions. Contrary to the standard implicit 
method by Fincham [Mol. Simul., 8, 165 (1992)], the revised angular 
momentum approach presented corresponds completely to the leapfrog 
idea on interpolation of dynamical variables without using any extrap- 
olations. The proposed scheme intrinsically preserves rigid molecular 
structures and considerably improves stability properties and energy 
conservation. As is demonstrated on the basis of simulations for wa- 
ter, it allows to reproduce correct results with extra large step sizes of 
order 5 fs and 10 fs in the cases of energy- and temperature-conserving 
dynamics, respectively. We show also that iterative solutions can be 
avoided within our implicit scheme shifting from quaternions to the 
entire rotation-matrix representation. 

Keywords: Numerical algorithms; Long-term integration; Motion of rigid bod- 
; Polyatomic molecules 



I. INTRODUCTION 



Computer experiment by the method of molecular dynamics (MD) is intensively 
exploited in solving various tasks of chemical physics [1], biochemistry [2] and bi- 
ology [3]. Among these are investigations of structure and dynamical properties 
of molecular liquids which normally are treated as collections of rigid bodies. De- 
spite the long prehistory of MD simulation, the development of efficient and stable 
algorithms for the integration of motion for such systems still remains an actual 
problem. 

Usually molecular movements are simulated using constrained dynamics [4-7] 
in which the phase trajectory of each atom is evaluated by Newton's equations, 
while the molecular structures are maintained by holonomic constraints to keep 
intramolecular bond distances. Although the atomic-constraint technique can be 
applied, in principle, to arbitrary polyatomics regardless of its chemical structure 
and size, it appears to be very sophisticated to implement for some particular models. 
For example, when there are more than two, three or four interaction sites per 
molecule for linear, planar and three-dimensional bodies, bond lengths and angles 
cannot be fixed uniquely [5]. Systems of point molecules with embedded multipoles 
present additional complexities too, because then the intermolecular forces cannot 
easily be decomposed into direct site-site interactions. The limitation of constrained 
dynamics is also caused by the fact that constraint forces are calculated at each time 
step of the produced trajectory to balance all other potential forces in the system. 
As the number of atoms in each molecule is increased, the number of constraints 
raises dramatically, resulting in a decreased speed of computations. Moreover, to 
reproduce the rigid molecular structure, cumbersome systems of nonlinear equations 
must be solved iteratively. This can lead to a problem for molecules with light 
hydrogen atoms or with linear or planar fragments. In this case, the algorithm 
converges rather slowly [8] already at relative small step sizes and, thus, it requires 
a considerably portion of the computational time. Recently, it was shown that a 
non- iterative calculation of constraint forces is possible [9] , but this is practical only 
for simple models in which the problem can be reduced to inversion of a banded 
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matrix [10, 11]. 

Some of the limitations just mentioned are absent in the molecular approach, 
when the displacements of rigid bodies are analyzed in view of translational and 
rotational motions. The translational dynamics is defined by motion of molecular 
centres of masses, whereas the orientations typically are expressed in terms of quater- 
nions [12-14] or principal-axis vectors [13]. The straightforward parameterization 
of orientational degrees of freedom, Euler angles, is very inefficient for numerical 
calculations because of singularities inherent in the description [12, 15, 16]. Mul- 
tistep predictor-corrector methods were applied to integrate rotational motion in 
early investigations [17-20]. As was soon established, the extra order obtained in 
these methods is not relevant, because the forces existing in a real system are not 
sufficiently smooth. As a result, high-order schemes appear to be less accurate at 
normal step sizes than low-order integrators, such as Verlet [21], velocity Verlet [22] 
and leapfrog [23] ones. The last algorithms are also the most efficient in view of 
cost measured in terms of force evaluations. That is why, they are widely used in 
different approaches, for instance, in the atomic-constraint technique, to integrate 
translational motion. These traditional algorithms were derived, however, assum- 
ing that velocities and forces are coordinate- and velocity-independent, respectively. 
In general, the time derivatives of orientational positions may depend not only on 
angular velocities but also on these positions themselves resulting in the explicit 
velocity-dependence of angular accelerations. Therefore, additional revisions are 
necessary to apply the standard integrators to rotational motion. 

In the atomic approach, the problem with the coordinate and velocity depen- 
dencies is circumvented by involving fundamental variables, namely, the individual 
Cartesian coordinates of atomic sites. Similarly, this problem can be solved within 
the molecular approach choosing appropriate generalized coordinates in orientational 
space. Ahlrichs and Brode proposed a method [24] in which the principal axes of 
molecules are treated as pseudo particles and constraint forces are introduced to 
maintain their orthonormality. Kol et al. considered the entire rotation matrix and 
the corresponding conjugate momentum as dynamical variables [25]. The rotation 
matrices can be evaluated within the usual Verlet or leapfrog frameworks, using 
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either recursive [24] or iterative [25] procedures, respectively. The recursive method 
behaves relatively poor with respect to long-term stability of energy, whereas the 
iterative procedure requires again, as in the case of constrained dynamics, to find 
solutions for systems of highly nonlinear equations. In general, the convergence of 
iterations is not guaranteed and looping becomes possible even at not very large 
step times. Examples for not so well behaved cases are models with almost linear or 
planar molecules, when the diagonal mass matrices are hard to numerical inversion 
since they have one or two elements which are very close to zero. The extension of 
the atomic and pseudo-particle approaches to temperature-conserving dynamics is 
also a difficult problem, given that the rigid-reactions and temperature-constraint 
forces are coupled between themselves in a very complicated manner. 

A viable alternative to integrate the rigid-body motion has been provided by 
explicit and implicit angular-momentum algorithms of Fincham [26-27]. This was 
the first attempt to adopt the leapfrog framework to rotational motion in its purely 
classical treatment. The chief advantage of these rotational leapfrog algorithms is 
the possibility to perform thermostatted simulations. However, even in the case of 
a more stable implicit algorithm, the total energy fluctuations in energy-conserving 
simulations are too big with respect to those identified in the atomic-constraint 
technique. Moreover, despite the fact that no constraint forces are necessary in 
the rotation dynamics, the rigidness of molecules is not satisfied automatically, be- 
cause the equations of motion are not solved exactly. Usually, the artificial rescaling 
method [19, 27] is used to preserve the unit norm of quaternions and, as a conse- 
quence, to ensure the molecular rigidity. Recently [28], it has been shown that the 
crude renormalization can be replaced by a more rigorous procedure introducing 
so-called numerical constraints. As a result, quaternion [28] and principal-axis [29] 
algorithms were devised within the velocity Verlet framework. It was demonstrated 
[29, 30] that these algorithms conserve the total energy better than the implicit 
leapfrog integrator [27], but worse with respect to the atomic-constraint method, 
especially in the case of long-duration simulations with large step sizes. 

Quite recently, to improve the stability, a new angular- velocity leapfrog algorithm 
for rigid-body simulations has been introduced [30]. The automatic preservation of 
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rigid structures and good stability properties can be related to its main advan- 
tages. But a common drawback, existing in all long-term stable integrators on rigid 
polyatomics, still remained here, namely, the necessity to solve by iteration the sys- 
tems of nonlinear equations. Although such equations are much simpler than those 
arising in the atomic and pseudo-particle approaches, the iterative solution should 
be considered as a negative feature. Moreover, since the nonlinear equations are 
with respect to velocities, it is not so simple matter to extend the angular-velocity 
algorithm to a thermostatted version. 

This study presents a modified formulation of the angular- momentum approach 
within the leapfrog framework. Unlike the standard approach by Fincham [27], the 
new formulation is based on more natural interpolations of dynamical variables and 
it uses no extrapolation. The algorithm derived appears to be free of all the draw- 
backs inherent in previous descriptions. It can easily be implemented to arbitrary 
rigid bodies and applied to temperature-conserving dynamics. The integrator ex- 
hibits an excellent energy conservation, intrinsically reproduces rigid structures and 
allows to avoid any iterative procedures at all. 



II. BASIC EQUATIONS OF MOTION 

Let us consider a system of N interacting rigid bodies. According to the classical 
approach, any movements of a body can be presented as the sum of two motions, 
namely, a translational displacement of the centre of mass and a rotation about 
this centre. The translational displacements in the system are expressed in terms 
of the centre-of-mass positions r-j and velocities Vj, where i = 1, . . . , N, given in a 
space-fixed laboratory frame. The time evolution of such quantities can described by 
writing Newton's law in the form of two per particle three-dimensional differential 
equations of first order, 

dvi 
di-j 

d^ = V - 
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where fj is the total force acting on body % due to the interactions with all the rest 
of particles and rrii denotes the mass of the body. 



A. Different forms of the equations for rotational motion 



To determine the rotational motion, one needs to use frames attached to each 
body together with the laboratory system of coordinates. It is more convenient for 
further consideration to direct the body-fixed-frame axes along the principal axes of 
the particle, which pass through its centre of mass. Then the matrix Jj of moments 
of inertia will be diagonal and time-independent in the body-fixed frame. We will use 
the convention that small letters stand for the representation of variables in the fixed 
laboratory frame, whereas their counterparts in the body frame will be designated 
by capital letters. The transitions E = A^e and e = A~ X E between these both 
representations of vectors e and E in the laboratory and body frames, respectively, 
can be defined by the 3x3 time-dependent rotation matrix Aj(t). Such a matrix 
must satisfy the orthogonormality condition Af Aj = I = AjA^~, or in other words 
A" 1 = Af, to ensure the invariance E + E = e + e of quadratic norms for vectors e 
and E. In our notations A -1 and A + are the matrices inversed and transposed to 
A, correspondingly, and I denotes the unit matrix. 

Let Aj be an arbitrary vector fixed in the body. By definition, such a vector 
does not change in time in the body-fixed frame, dAj/dt = 0. The angular velocity 
u?j is introduced differentiating its counterpart Si(t) = A+(i)A, in the laboratory 
frame over time, dSi/dt = u>iX<5j. Then, using the equality aj>jX<5j = W + (a;j)<5j 
and the orthonormality of Aj, the rate of change in time of the orientational matrix 
can be expressed in terms of either laboratory u>j or principal fij = AjO>j angular 
velocity as 




dA, 



(2) 



where 
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is a skewsymmetric matrix, i.e., W + (f2;) = — W(f2j), and Q x , Vt % Y and Vt l z are 
components of vector f2j. 

From the orthogonormality condition it follows that maximum three indepen- 
dent parameters are really necessary to describe orientations of a rigid body and to 
evaluate the nine elements of the rotation matrix. However, the well-known param- 
eterization of A, in terms of three Eulerian angles [12] is unsuitable for numerical 
calculations because of the singularities. In the body- vector representation [13, 24, 
25, 29], all the elements of the rotation matrix A« are considered as dynamical vari- 
ables. These variables present, in fact, Cartesian coordinates of three principal axes 
XYZ of the body in the laboratory frame. The alternative approach applies the 
quaternion parameterization [4, 9] of rotation matrices, 



A(qi 
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-2(&7< + QXi) & - Vi ~ C! + xl 2(ViXi ~ M 

livid-tiiXi) -2(tii(i + ViXi) -g-Vi+Ci+x 2 

where = rji, Q, Xi) + is a vector-column consisting of four quaternion compo- 
nents. Using the normalization condition q^~qi = £f + rjf + (f + x! = 1) which 
ensures the orthonormality of Aj(t) = A[qj(i)], the time derivatives of quaternions 
can be cast [13, 27, 30] in the form 
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y * L z 
where Q(ftj) is a skewsymmetric matrix again. 

Expressions (2) and (5) are rotation-motion analogues of the second line of 
Eq. (1) in the case of body-vector and quaternion representations, respectively. 
They must be complemented by equations defining the time evolution of angular 
velocities. The simplest form of these equation is obtained for the angular momenta 
lj = A.f'Li of bodies in the laboratory frame, where Lj = Jjfij are principal angular 
momenta. The result is 



where kj is the torque exerted on body i with respect to its centre of mass. The angu- 
lar velocities in the body- or space-fixed frames can easily be reproduced, whenever 
they are needed, using the relations f2j = J~ 1 A i lj and u>j = A+ f2j = jj rl L, where 
ji = AfJiAi is the time-dependent matrix of moments of inertia in the labora- 
tory frame. Another way lies in involving explicit equations for principal angular 
velocities. Such equations, known also as Euler's ones, can be derived substitut- 
ing lj = AfJifli into Eq. (6) and using equations of motion (2) for orientational 
matrices. As a result, one obtains 



where Kj = Ajkj are the principal torques. Formally replacing the quantities fij, Kj 
and Ji by Wj, kj and jj yields quite similar equations of motion for angular velocities 
u>i in the laboratory frame. 

It is worth remarking that the body- vector (Eq. (2)) and quaternion (Eq. (5)) 
representations as well as the angular-momentum (Eq. (6)) and angular- velocity 
(Eq. (7)) approaches are completely equivalent between themselves from the math- 
ematical point of view. For numerical evaluations, the preference must be given to 
equations which allow to be integrated in the simplest manner with the greatest 
precision and the best stability. In the present study we shall deal with more simple 
equations of motion (6) for angular momenta in the laboratory frame rather than 
with equations (7) for principal angular velocities. In such a way, difficulties with 
the velocity-dependence of angular accelerations are excluded automatically. More- 
over, we shall show that the angular-momentum approach allows to obviate iterative 
solutions within a leapfrog framework choosing the entire-rotation-matrix elements, 
instead of quaternions, as orientational variables. Thus, the body-vector represen- 
tation should be considered as a more preferable method for such calculations. 

III. THE REVISED ANGULAR-MOMENTUM APPROACH 

In the case of translational motion, equations (1) can readily be integrated with 
the help of the usual [23] leapfrog algorithm: 




(7) 
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Vi(t +*)=v i (t-$)+h fi(t)/m + 0(h 3 ) , 
Ti(t + h)= Ti {t) + hvi(t + |) + 0{h 3 ) , 



(8) 



where h is the time increment, and forces fj(t) are evaluated using known spatial 
coordinates Ti(t) and A«(t). The truncation local errors, appearing during such an 
integration, are of order h 3 in both coordinates and velocities. If an estimator of 
Vj(t) is required, for example to evaluate the total energy, the usual choice is 



where interpolation uncertainties 0(h 2 ) are in the self-consistency with the second 
order of global errors (one order lower than that for local errors) of the leapfrog 
integrator (8). 



For the rotational motion the time derivatives of orientational positions (Eqs. (2) 
and (5)) depend not only on angular velocities but also on these positions them- 
selves. This difficulty cannot be handled with a simple leapfrog scheme in which the 
positions and velocities are known at different times. Relatively recently, Fincham 
[27] has proposed a solution to the problem by introducing an implicit leapfrog-like 
algorithm. His method can briefly be described as follows. 

First, quite analogously to the case of translational- velocity evaluations (fist line 
of Eq. (8)), angular- momentum equations (6) are integrated as 



At this stage the principal angular velocities fli(t) can be calculated using the rela- 
tion 



v i (t) = l[v i (t-§)+v i (t+§)]+0(h 2 ), 



(9) 



A. Standard rotational leapfrog algorithm 



M* + f) = h(t-!i) + hk i {t) + o(h 3 ). 



(10) 




m(t) 



(ii) 



and the propagation 



W) = I [W - 1) + h(t + 1)] = h(t - 1) + 



(12) 
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of angular momenta to on-step level of time. 

Further, according to the leapfrog framework, the evaluation of orientational 
coordinates must be performed as 

Si(t + h) = Si(t) + hSi(t + |) + 0(h 3 ) , (13) 

where Sj = dSj/dt = H(f2j)Sj, and either S« = Aj and H = W or Sj = q 8 and 
H = Q in the case of either entire-matrix or quaternion space, respectively. Note 
that in the quaternion representation the orientational matrices Aj(i) = Ajfq^t)] 
appear implicitly, and they are computed via relation (4) using quaternion values 
q,i(t). As far as the quantities Sj and f2j are not known at mid-step level t + |, 
it was assumed to propagate the time derivatives of S« by means of the relation 
Si(t + |) = |[Si(f) + Si(t + h)} + 0(h 2 ), i.e., 

Si(t + §) = ± [H(fi i (*))S i (*) + H(a(t + h))Si(t + h)] + 0(h 2 ) , (14) 

where 

fii(t + h) = J^Mt + h)\(t + h) . (15) 

Propagation (14) requires in its turn the knowledge of advanced angular momenta 
li(t + h) which were predicted by writing 

l(t + h) = \{t + |) + |k,(0 + 0(h 2 ) . (16) 

In view of (14) and (15), relation (13) is an implicit system of equations with 
respect to elements of Si(t + h), defined through the auxiliary parameters lj(t) and 
li(t + h) which are not stored, but used to calculate the angular velocities fti(t) 
and Qi(t + h) in the body frame. The system can be solved by iteration taking 
Sf\t + h) = St(t) + hU(ni(t))Si(t) as the initial guess. 

A thermostatted version is based on interpolations (9) and (12) of on-step trans- 
lational velocities and angular momenta. Such interpolations are used in micro- 
canonical simulations to evaluate the kinetic temperature T(t) = T({vj(t), fij(i)}) = 
^Eiljmv^) +Ea X ' Z JLK 2 (t)}, where J l xx , J l YY and .P zz are nonzero ele- 
ments of matrix Jj, /c B is the Boltzmann's constant and I = 6 denotes the number 
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of degrees of freedom per particle (for linear bodies 1 = 5). This allows to syn- 
chronize in time the temperature with the potential energy U(t) = £/({r;(t), Sj(t)}) 
and, therefore, to calculate the total energy E(t) = l -^ s -T(t) + U(t) of the sys- 
tem. In the temperature-conserving dynamics, on-step velocities and angular mo- 
menta are modified as v-(t) = /3(t)vi(t) and Y^t) = (3(t)li(t) using the scaling factor 
f3(t) = y / T /T(t), where T is the required constant temperature [27, 31]. The 
velocity integration is completed by 

v:(t + |) = [2-r 1 (t)K(t) + |f i (t)/m, (17) 

Y t (t+D = [2-p- i (t)]m + ik t (t) (is) 

which satisfy the interpolations v-(£) = |[v;(t - §) + v-(i + §)], l-(t) = - 
|) + l^(t + |)] and the constant-temperature condition T({v^(t), f2J(i)}) = T , where 
fi-(t) = J l ~ 1 Aj(t)l-(t). Finally, the translational and orientational positions are 
updated according to the same equations replacing v;(t + |) by v-(t + |) and lj(t+|) 
by + 

B. Revised leapfrog algorithm 

As has been established [29, 30], the rotational leapfrog algorithm, described in 
the preceding subsection, exhibits rather poor long-term stability of energy with 
respect to atomic-constraint integrators [4-7] , for example. Moreover, it requires it- 
erative solutions and does not conserve the unit norm and orthonormality of quater- 
nions and orientational matrices. For this reason, a question arises how about the 
existence of a revised scheme which is free of all these drawbacks and which has all 
advantages of the standard approach. We shall show now that such a scheme really 
exists. 

First of all, one points out some factors which can explain bad stability properties 
of the standard scheme. When calculating orientational variables, the Fincham's 
algorithm uses up three additional estimators, namely, the propagations for on-step 
angular momentum l^t) (Eq. (12)) and mid-step time derivative Sj(t + |) (Eq. (14)) 
as well as the prediction (Eq. (16)) of angular momentum l^t+h). Among these only 
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the first two evaluations can be classified as interpolations which correspond to a 
simple averaging over the two nearest neighbouring values. At the same time, the last 
prediction (16) presents, in fact, an extrapolation that is, strictly speaking, beyond 
the leapfrog framework. Indeed, applying equation (12) for the next step time 
t = t + h yields the following interpolated values li(t + h) = L,(t + |) + |kj(t + h) for 
angular momenta, which differ from previously predicted ones, i.e., li(t+h) ^ \i{t+h) 
and, as a consequence, Cli(t + h) ^ fl^t + h). Extrapolations are commonly used in 
low-precision explicit schemes and they should be absent in more accurate implicit 
integrators. 

The main idea of the revised approach is to derive an implicit equation for 
Si(t+h) reducing the number of auxiliary interpolations to a minimum and involving 
no extrapolations. This can be realized starting from the same evaluation (10) for 
mid-step angular momenta, but treating the time derivatives Sj(t + |) = H(f2j(£ + 
|))Si(t + |) in a somewhat other way. As was mentioned earlier, these derivatives 
are necessary to evaluate orientational positions (Eq. (13)), and they require the 
knowledge of two per body quantities, namely, + |) and S { (t + |). It is crucial 
to remark that since the advanced angular momenta h(t + |) are already known, 
such two quantities are not independent but connected between themselves by the 
relation 

a(t+|)=J- 1 A i (t + |)l l (t + |). (19) 

Then, as can be seen easily, the calculation of Sj(t+ 1) is reduced to a propagation of 
one variable only, namely, Sj(t + 1). It is quite naturally to perform this propagation 
by writing 

S»(* + |) = \ [Si(t) + Si(f + h)] + 0{h 2 ) (20) 
and the algorithm proceeds as follows 

Si(t + h) = Si(t) + hu(n i (t + f )) Si(t + 1) + o(h 3 ) . (21) 

Taking into account expressions (19) and (20), matrix equation (21) constitutes 
an implicit system for unknown elements of Sj(t + h). As for the usual scheme, 
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the system can be solved iteratively, putting initially S- '(t + h) for Sj(t + h) in all 
nonlinear terms collected in the right-hand side of (21). Then the obtained values for 
Si(t + h) in the left-hand side are considered as initial guesses for the next iteration. 
The convergence of iterations is justified by the smallness of nonlinear terms which 
are proportional to the step size h. 

In such a way, we have derived a new leapfrog algorithm to integrate orienta- 
tional degrees of freedom. It involves only one auxiliary interpolation (20) which 
is completely in the spirit of the leapfrog framework. Moreover, this interpolation 
concerns the most slow variables Sj, rather than their more fast time derivatives Sj 
and angular momenta lj, thus, leading to an increased precision of the calculations. 
When on-step temperature T(t) is required, for instance to check the energy conser- 
vation, we can apply usual interpolation (12) of angular momenta and relation (11) 
for velocities. It is worth underlining that, unlike the standard rotational integrator, 
the angular-momentum interpolation errors are not introduced into trajectories (21) 
produced by the revised algorithm at least within the energy-conserving dynamics. 

The extension of the revised scheme to a thermostatted version is trivial. Using 
the calculated temperature T(t) we define the scaling factor (3(t) = \JT /T(i). The 
mid-step angular momenta lj(i+|) are then replaced by their modified values l-(t+|) 
(see Eq. (18)) and substituted into Eq. (19) to continue the integration process 
according to equations (20) and (21). 

Besides the evident simplicity of the revised approach with respect to the stan- 
dard scheme, a very nice surprise is that the unit norm of quaternions and the 
orthonormality of orientational matrices appear to be now by numerical integrals of 
motion. Indeed, considering the quantity Cli(t + |) as a parameter and explicitly 
using coordinate interpolation (20), we can present Eq. (21) in the equivalent form 

Si(t + h) = [I - |H,(t + l)]- 1 ^ + |H,(t + !)]Si(t) , (22) 

where Hj(i + |) = H(f2j(i + |)) and it is understood that I designates either three- 
or four-dimensional unit matrix in the principal-axis or quaternion domain, respec- 
tively. It can be checked readily that the matrix (I — 0) _1 (I + 0) is orthonormal for 
an arbitrary skewsymmetric matrix @ + = — 0. As far as the matrix H is skewsym- 
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metrical by definition, the following important statement emerges immediately. If 
initially the orthonormality of Sj(£) is fulfilled, it will be satisfied perfectly for the 
advanced matrices Sj(£ + h) as well, despite an approximate character of the inte- 
gration process. Thus, no artificial or constraint normalizations and no recursive 
procedures are necessary to conserve the rigidness of molecules. 

The alternative presentation (22) may be more useful for iterating since it pro- 
vides the orthonormality of Si(t + h) at each iteration step and leads to an increased 
speed of the convergence. Because of this, we show Eq. (22) more explicitly, 

*(* + Q = 1[l '^it^t kCtl *(*) - G,(t, h) , (23) 

h 2 r,2(+ i h\\ i u\\T i h' 



I [1 - ^(t + |)] + hW t + \v t 



Mt + h) = 4 ; V * ' 2 Mt) = B t (t, h) AM , (24) 



for the cases of quaternion and entire-rotation-matrix representations, respectively, 
where expressions (3) and (5) for matrices W< = W(fJ;(t + §)) and Q; = Q(Qi(t + 
|)) have been taken into account, Gi(t,h) and T>i(t,h) are orthonormal evolution 
matrices, and [Pj] Q| g = denotes a symmetric matrix which, like Wj and Qj, 

is calculated using principal angular velocities (19). In view of the equalities Wf = 
Pj — Ofl and Qf = — the evolution matrices can be cast also in the matrix- 
exponential forms 



Gi(t, h) = exp[0jQj/r2j] t+ h , fa = 2arcsin 



fq(t + f) 

i + gn?(* + J) 



Di(t, /i) = exp[(/?jWj/Qi] t+ fe , pi = arcsin J — - - 



(25) 



i + m 2 (t+|) 



Then it becomes clear that the matrices Dj and Gj define three- and four- 
dimensional rotations on angles ifii and fa in the laboratory frame and quaternion 
space, respectively. In the first case the rotation is carried out around the unit vec- 
tor fli/Qi\ t+ h, whereas in the second one it is performed around an orth which is 
perpendicular to all four orths of quaternion space. 
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C. Avoidance of iterative solutions 



Another excellent feature of the algorithm is that within the entire matrix rep- 
resentation, equation (24) can be handled in a non-iterative way using so-called 
quasianalytical solutions for mid-step angular velocities fti(t + |). To show this we 
first perform a set of further transformations. Remembering that now Sj = Aj and 
Hj = Wj, one adds the matrix Aj(t) to the both sides of Eq. (21) and divides the 
obtained equation by factor 2. Then using coordinate propagation (20) leads to 

Mt + |) = A,(t) + |W(^(t + I)) Mt + |) . (26) 

Multiplying Eq. (26) on the matrix J" 1 from the left and additionally on the vector 
lj(t + |) from the right, and taking into account definition (19) yields 

n t (t + 1) = jr'MtMt + §) + fj-iw^t + f))JA(* + |) • (27) 

Therefore, the iterative problem is much simplified, because it is reduced to 
finding solutions to three-dimensional vector equation (27) for three unknown com- 
ponents fix, fiy and Vl z of f2;(t + |) rather than to matrix equation (24) (or (21)) 
for nine unknowns elements of Aj(t + h). Equation (27) can be solved iteratively 
again, choosing J^ 1 Aj(t)lj(t + |) as the initial guess for fli(t + |). 

A next simplification lies in the following. Let us rewrite equation (27) in the 
explicit form 

fix = 9x + hg x n Y ^lz , 

tt Y = 9 Y + hQ Y VL z VLx , (28) 
fiz = 0z + IiqzQxQy , 

where g x = (J YY - 4z)/( 2 4x)» Qy = (4z ~ Jxx)/( 2 Jyy), Qz = {Ax ~ 
Jyy)/{2Jzz) = ~(@x + Qy), and 6 x ,y,z are the components of known vector 
J~ 1 Aj(t)lj(t+ |), keeping in mind that vector li(t+ |) must be replaced by l-(t+ |) 
in the case of temperature-conserving dynamics. Unless J l xx ^ J YY ^ J l zz , the 
system of equations (28) appears to be linear and, therefore, it can easily be solved 
exactly (see subsect. III. D, where specific models are described). Here, we consider 
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the most general case when all the principal moments of molecules are different and 
assume for definiteness that J xx < J YY < 3 % zz . Then the first two unknowns £l x 
and f2y are the most fast variables and they should be excluded from the iteration 
to increase the convergence. Such an excluding indeed can be realized solving the 
first two equations of (28) with respect to Q x and Q Y . The result is 

e x + hQx8 Y n z 6 Y + hg Y 9 x n z , on , 

Qx ~ i + h 2 v 2 n z 2 ' Qy ~ i + h 2 v 2 n z 2 ' (29) 

where < v 2 = —QxQy < 1/4. The last inequalities follow from the requirements 
Jaa > and J aa < Jpp + J 77 imposed on principal moments of inertia, where 
(a, /5, 7) denote an array of three cyclic permutations of (X,Y,Z). In view of (29), 
only the third equation of system (28) really needs to be iterated with respect to 
one variable Q z . Since Q z is the most slow quantity, a well convergence can be 
guaranteed even for not so well normally behaved case as an almost linear body, 
when J l xx <C J YY < J zz - 

Finally, one considers the question of how to obviate iterative solutions at all. 
Substituting the result (29) into the third equation of system (28) and presenting the 
Zth component of the angular velocity in the form fl z = s + 5 yields the following 
algebraic equation 

a + ai8 + a 2 5 2 + a 3 5 3 + a A 5 A + a 5 5 5 = (30) 

with the coefficients 

a = (so - 6 z )$l - hQ Z [0x0 Y fi- + h(g Y 6 x + g x Y )s o } , 

cn = #+ - h 2 {(g Y 6 x + g x e 2 Y )g z - u 2 s [(5s - M z )d + + 1hQ x Q Y g z \} , 

a 2 = h 2 u 2 [Qs - 26 z + hg z 6 x 6 Y + h 2 u 2 s 2 (10s - Q6 Z )} , (31) 

a 3 = 2h 2 u 2 [l + h 2 u 2 s (5s - 26 z )] , 

04 = /i 4 z/ 4 (5so _ 6z) , Q<5 = h 4 is 4 , 

where ~&± — 1 ± }{ 2 v 2 s\. The equation (30) is fifth order and the corresponding solu- 
tions for Q z are independent on parameter so, provided the unknown 5 is precisely 
determined. However, as is well known, only algebraic equations of fourth or less 
order allow to be solved in quadratures. 
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The main idea of quasianalytical solutions consists in the fact that actual MD 
simulations are performed with relative small values of the time step h. Then it is 
necessary to choose the parameter s as a good prediction for Q z to be entitled to 
ignore high-order terms in the left-hand side of Eq. (30). The simplest choice for 
this can be found putting Qx^y = Qx®y + 0(h) in the right-hand side of the third 
equation of system (28). As a result, one obtains 

s = Z + h Qz 6 x e Y (32) 

that represents the original value of Q z with second-order truncation errors, so that 
5 = 0(h 2 ). It is easy to see that in this case the two last terms a 4 <5 4 and a 5 5 5 in the 
left-hand side of Eq. (30) behaves as 0(h 12 ) and 0(h u ), respectively. Taking into 
account the smallness of h, such terms can merely be omitted without any loss of the 
precision, because they involve uncertainties of order 0(h 12 ) into the solutions and 
appear to be too small with respect to third-order truncation errors 0(h 3 ) inherent 
initially in the algorithm. 

Eq. (30) is now transformed to the third-order algebraic equation 

a + axS + a 2 5 2 + a 3 5 3 = 0(h 12 ) (33) 

which can easily be solved analytically. The result is 

S 1 = -la 2 /a 3 + c-b/c + 0(h 12 ), (34) 
£2,3 = -|a 2 /a 3 -\[c-b/c± iV3(c + b/c)\ + 0{h 12 ) , 

where 

b = |(3aia 3 - a 2 2 )/a 2 3 , 

c=(p+^/b 3 +p 2 f 3 , (35) 
p = ^(9aia 2 a 3 - 27a aj - 2a 3 2 )/a 3 3 . 

Among three solutions (34), only the first one 5\ is real and satisfies the physical 
boundary condition ~ h 2 when h goes to zero (the other two solutions #2,3 are purely 
imaginary at h — > and they tend to infinity as ~ ±i/h). 

Therefore, the desired Zth component of the angular velocity is 
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&>z = so + Si . 



(36) 



The rest two components fix and Qy are reproduced on the basis of equalities (29). 
The obtained vector fl^t + |) = (fyx-,^V,^z) is substituted into equation (24) to 
perform the explicit evaluation for advanced orientational matrices A^t + h). This 
completes the algorithm. 

D. Implementations for particular models 

There are two main classes of models for interacting rigid bodies, which most 
frequently are applied in MD simulations. More realistic so-called interaction site 
models can be related to the first class. For these models, each ith body of the system 
is considered as a molecule which in its turn is composed of Mj point interaction sites 
(atoms). The rigid structure of molecules is completely defined by time-independent 
vector-positions A" (a = 1, . . . , Mj) of atom a within molecule % in the body frame, 
whereas these positions in the laboratory frame are: r"(t) = r-j(t) + A+(t)A". Using 
the known site-site potentials the desired molecular forces and torques can 

easily be computed as f, = E^3)2> ~ r 'D and k * = ^5)5>( r * " r xf S'> 

respectively, where = —dufj/drfj and r°* = r" — r b -. The second class is point 
molecules (maxj ;a)fe |A" — A£| — >• 0) with embedded multipoles. The most popular 
model belonging this class is a system of point electro-dipoles. The molecular forces 
and torques caused by dipole-dipole interactions can be calculated using the relations 

fi = Efij^i) £[ r v{AVA*j - ;?:(/V r ij)(/V r v)} + /x^Aij-Tij) + fi> ! ' r «i)] and k; = 
7r/*,x;-4r, ; (// ;; T; ; ) - fij], respectively, where = - and ^ denotes 
the dipole moment of ith molecule. 

Although the proposed algorithm can be implemented for arbitrary rigid mod- 
els, some simplifications with respect to the general formulation are possible us- 
ing special properties of the body. The simplest case is bodies with a spherical 
distribution of mass, when all the moments of inertia are equal between them- 
selves, i.e., when J l xx = Jy Y = J zz = Jj and, thus, Jj = JT = jj. Then it 
is more convenient to work with equations (2) for rotational matrices, presented 
in terms of angular velocities cjj = j" 1 ^ = h/Ji in the laboratory frame, i.e., 
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with dAj/dt = AjW(uJi). The leapfrog trajectories for these equations are ob- 
vious: Ai(t + h) = A^^exp^W^Wj)/^]^!, where u>i(t + §) = h(t + \)/Ji and 
(fi = arcsin[/icjj/(l + ^uf)] t+ h. 

For some particular models, the orientational part of the intermolecular potential 
can be expressed using only unit three-component vectors p t passing through the 
centres of mass of molecules. The examples are point dipole interactions, when 
p i = pj or when all force sites of the molecule are aligned along p iy resulting in 
torques which are perpendicular to p iy i.e., k i -p i = 0. If then additionally J« = JJ 
(for the last example this can be possible when forceless mass sites are placed in such 
a way to ensure the spherical mass distribution), it is no longer necessary to deal 
with orientational matrices or quaternions. In this case the equation for p i looks as 
dpj/dt = W + (u> i )p i with the solution p^t + h) = exp[— (piW^aj^/ui^+hp^t). 

For molecules with cylindric distribution of mass sites, when two of three of prin- 
cipal moments of inertia are equal, the numerical trajectory can also be determined 
in a simpler manner. Let us assume for defmiteness that J xx = J YY ^ J l zz and 
J zz 7^ 0. Then arbitrary two perpendicular between themselves axes, lying in the 
plane perpendicular to the Zth. principal axis, can be considered initially as X- and 
Y-th. principal orths. Since now q z = 0, the Zth component of the angular velocity 
is found automatically, Qz — Qz- As in the general case, the two rest solutions fix 
and fl Y of system (28) are calculated on the basis of Eq. (29) taking into account 
that Qy = —Qx, whereas the orientational matrices are evaluated via Eq. (24). 

A special attention must be paid for purely linear molecules when J xx = J YY 7^ 
J zz = and each body has two, instead of free, orientational degrees of freedom. 
The relative positions r"(t) — rj(t) = A^p^t) of all atoms within a linear molecule 
can be expressed in terms of an unit vector p { and besides ki-p { = one finds that 
\Lj\ z = J zz fl l z = 0. The rotational part |( J XX VT X 2 + Jy y f2y 2 ) of the kinetic energy 
is also indifferent to the Zth. component Q z of the principal angular velocity. Such 
a component causes irrelevant rotations of the molecule around p r axis and it does 
not lead to any change of and the potential energy. It can be shown that the 
angular-momentum approach allows to reproduce the correct time evolution of two- 
dimensional unit vector p i by the three-dimensional leapfrog rotation p^t + h) — 
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exp[(/?jWi(f2j)/Jli] t+ |Pj(t) putting formally VL l z = 0, where two other components 
of f2«(t + §) are fix = X and f2y = 9 Y (this immediately follows from Eq. (28)). 
Planar molecules do not present a specific case within our approach and they are 
handled in the usual way as tree-dimensional bodies. 

IV. NUMERICAL VERIFICATION OF THE ALGORITHM 

The system chosen for numerical tests was the TIP4P model (M = 4) of water 
[32] at a density of mN/V = 1 g cm" 3 and a temperature of T = 298 K. Because of 
the low moments of inertia of the water molecule and the large torques due to the 
site-site interactions, such a system should provide a very severe test for rotational 
algorithms. In order to reduce cut-off effects to a minimum we have applied an inter- 
action site reaction field geometry [33] and a cubic sample with N = 256 molecules. 
All runs were started from an identical well equilibrated configuration. The MD sim- 
ulations have been carried out in both energy- conserving (NVE) and thermostatted 
(NVT) ensembles. The equations of rotational motion were integrated using the 
standard quaternion integrator [27] and our revised leapfrog algorithm. As far as 
water is usually [34] simulated in an NVE ensemble by the atomic-constraint tech- 
nique [4, 5], the corresponding calculations on this approach and the angular- velocity 
Verlet method [28] were performed for the purpose of comparison as well. All the 
approaches required almost the same computer time per step given that near 97% 
of the total time were spent to evaluate pair interactions. 

The following thermodynamic quantities were evaluated: total energy, potential 
energy, temperature, specific heat at constant volume, and mean-square forces and 
torques. The structure of the TIP4P water was studied by determining the oxygen- 
oxygen and hydrogen-hydrogen radial distribution functions (RDFs). Orientational 
relaxation was investigated by evaluating the molecular dipole-axis autocorrelations. 
Centre-of-mass and angular-velocity time autocorrelation functions were also found. 
To reduce statistical noise, the measurements were averaged over 20 000 time steps. 

In the case of NVE dynamics to verify whether the phase trajectories are pro- 
duced properly, we applied the most important test on conservation of total energy 
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E of the system. The total energy fluctuations £ = [((E - (E)) 2 )] 1 / 2 /\(E)\ as func- 
tions of the length of the simulations over 10 000 time steps are plotted in Fig. 1 
(a)-(d) at four fixed step sizes, h — 1,2,3 and 4 fs. Both principal-axis (the boldest 
curves in subsets (a)-(d)) and quaternion representations were used to integrate the 
equations by the revised leapfrog algorithm. It has been established that the func- 
tions £ corresponding to these representations are practically the same. For this 
reason and to simplify the graph notations the results obtained within quaternion 
variables are shown (as crosses) only in subset (d) of the figure. 

As can be seen easily, the standard rotational leapfrog algorithm exhibits rel- 
atively bad stability properties and conserves the energy rather poor even at the 
smallest step size considered. It is worth remarking that investigating the system 
during shorter time periods with small step sizes, for example over 1000 time steps 
with h — 1 fs, one may come to very misleading conclusions on the energy con- 
servation. A significantly better pattern is observed for the angular velocity Verlet 
integrator. However, the improvements in stability are quite insufficient especially 
for moderate and large step sizes (h > 3 fs, subsets (c)-(d)). Finally, we can talk 
about the best energy conservation and long-term stability for the atomic-constraint 
scheme and the revised leapfrog algorithm which lead to virtually identical results. 

It is necessary to emphasize that within the principal-axis representation, the 
revised leapfrog trajectories were evaluated using the non-iterative quasianalytical 
scheme. The exact solutions (by means of iterations of Eq. (28)) were computed too 
to compare it with quasianalytical values. No deviation between both trajectories 
has been found up to h = 10 fs. They differed on each step by uncertainties of 
order round-off errors, so that the quasianalytical hypothesis appears to be in an 
excellent accord. At the same time, the quaternions converged at each step to a 
relative tolerance of 10~ 10 in average from 6 to 14 iterations with varying the step 
size from 2 fs to 6 fs. 

We also tried to avoid iterative procedures for quaternion variables by applying 
a hybrid leapfrog scheme when the quasianalytical solutions for mid-step angular 
velocities are substituted directly into orthogonormal matrices for quaternion eval- 
uation (23). However, the hybrid scheme leads to a significant loss of the precision 
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(see the long-dashed curve in subset (d) of Fig. 1). This is so because equation (27) 
for angular velocities was obtained on the basis of interpolation (20) for principal- 
axis variables, i.e., when A;(t + |) = |[Aj(£) + Aj(t + h)]. Using these veloci- 
ties in the quaternion space causes an inconsistency of such an interpolation with 
the corresponding interpolation q_i(t + |) = |[qi(i) + fk(t + h)] for quaternions 
since Ajq^i + |)] ^ |(Aj[q.j(t)] + Aj[qj(t + h)]). Therefore, to follow rigorously 
the leapfrog framework, the auxiliary mid-step angular velocities must be involved 
within the principal-axis representation exclusively. 

No shift of the total energy and temperature was observed during the revised 
leapfrog trajectories at h < 5 fs over a length of 20 000 time steps. As is well 
known, to reproduce features of an NVE ensemble correctly, the ratio T = S/U of 
the total energy fluctuations to the corresponding fluctuations U of the potential 
energy must be within a few per cent. The following levels of £ at the end of the 
revised leapfrog trajectories have been obtained: 0.0016, 0.0066, 0.017, 0.030, 0.051 
and 0.11 %. They correspond to T « 0.29, 1.2, 3.0, 5.4, 9.1 and 20 % at h = 1, 2, 3, 
4, 5 and 6 fs, respectively, where it was taken into account that U « 0.56% for the 
investigated system. The deviations in all the rest measured functions with respect 
to their benchmark values (obtained in the atomic-constraint NVE simulations with 
h — 2 fs) were in a complete agreement with the corresponding relative deviations 
T in the total energy conservation. For example, the results of the revised leapfrog 
algorithm at h — 2 fs were indistinguishable from the benchmark ones, whereas they 
differed as large as around 5%, 10% and 20% with increasing the time step to 4 fs, 
5 fs and 6 fs, respectively. However, the differences were much smaller than in the 
case of the standard rotational integrator. We see, therefore, that step sizes of order 
5 fs are still suitable for precise NVE calculations. Even a time step of 6 fs can be 
acceptable when a great precision is not so important, for instance, to achieve an 
equilibrium state. 

What about the NVT simulations? It is well established [27, 35] that ther- 
mostatted versions allow to perform reliable calculations with significantly greater 
step sizes than those used within the energy- conserving dynamics. To confirm such 
a statement, we have made NVT runs on the basis of our non-iterative revised 
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leapfrog algorithm (within principal-axis variables) and a thermostatted version of 
the standard implicit integrator [27] of Fincham. 

The oxygen-oxygen and hydrogen- hydrogen RDFs, calculated during the revised 
leapfrog trajectories for three different step sizes, h = 2, 8 and 10 fs, are plotted 
in Fig. 2a by the curves in comparison with the benchmark result (open circles). 
Note that the RDFs corresponding to h — 4 and 6 fs coincide completely with those 
for h = 2 fs and they are not included in the graph. A similar behaviour of RDFs 
was identified for the standard rotational integrator, but the results are somewhat 
worse especially at h — 8 and 10 fs. No drift of the potential energy was observed 
at h < 10 fs and h < 6 fs for the revised and standard algorithms, respectively. 
From the above, we can conclude that the revised leapfrog integrator is suitable for 
simulating even with huge step sizes of 10 fs, because then there is no detectable 
difference in RDFs. Other thermodynamic quantities such as the centre-of-mass and 
angular-velocity time autocorrelation functions appeared to be also close to genuine 
values. Quite recently, it was shown [36] that the time interval of 10 fs should 
be considered as an upper theoretical limit for the step size in MD simulations on 
water. We see, therefore, that this limit can be achieved in practice using the revised 
leapfrog algorithm. 

The molecular dipole-axis time autocorrelation function is the most sensitive 
quantity with respect to varying the step size. Such a function obtained within 
the standard (S) and revised (R) schemes at five fixed step sizes, h — 2, 4, 6, 8 
and 10 fs, is presented in Fig. 2b. For h < 6 fs the results of S- and R-schemes 
are indistinguishable between themselves. With increasing the step size to 8 fs or 
higher we can observe a systematic discrepancy which is smaller in the case of the 
R-scheme. Reliable results can be obtained here at time steps of h < 8 fs for both the 
standard and revised schemes. However within the standard approach, the solutions 
converged too slow already at h — 6 fs and they began to diverge at greater step sizes. 
To perform the simulations in this case, special time-consuming transformations to 
unsure the convergence have been applied. For the revised integrator which is free 
of iterations, the computer time did not depend on the step size. 
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V. CONCLUSIONS 



During last years there was a slow progress in the improvement of existing MD 
techniques concerning the numerical integration of motion for systems with inter- 
acting rigid bodies. We have attempted to remedy such a situation by formulating 
a revised angular- momentum approach within the leapfrog framework. As a result, 
a new integrating algorithm has been derived. The revised approach reduces the 
number of auxiliary interpolations to a minimum, applies the interpolations to the 
most slow variables and avoids any extrapolations. This has allowed to achieve the 
following two significant benefits: (i) all final expressions are evaluated explicitly 
without involving any iterative procedures, and (ii) the rigidity of bodies appears 
to be a numerical integral of motion. Another positive feature of the algorithm is 
its simplicity and universality for the implementation to arbitrary rigid structures 
with arbitrary types of interactions. 

As has been shown on the basis of actual simulations of water, the proposed algo- 
rithm exhibits very good stability properties and conserves the total energy in micro- 
canonical simulations with the same precision as the cumbersome atomic-constraint 
technique. In the case of temperature-conserving dynamics, reliable calculations are 
possible with huge step sizes around 10 fs. Such sizes are very close to the upper 
theoretical limit and unaccessible in usual approaches. 
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Figure captions 



Fig. 1. The total energy fluctuations as functions of the length of the NVE 
simulations on the TIP4P water, evaluated in various techniques at four fixed time 
steps: (a) 1 fs, (b) 2 fs, (c) 3 fs and (d) 4 fs. 

Fig. 2. Oxygen-oxygen (O-O) and hydrogen-hydrogen (H-H) radial distribution 
functions (a), and orientational relaxation (b), obtained in NVT simulations on the 
TIP4P water using the revised ((a), (b)) and standard ((b)) leapfrog algorithms. 
The results corresponding to the step sizes h — 2, 8 and 10 fs are plotted by bold 
solid, short-dashed and thin solid curves, respectively. Additional long-short dashed 
and dashed curves in (b) correspond to cases of h — 4 and 6 fs. The sets of 
curves related to standard and revised integrators are labelled in (b) as "S" and 
"R", respectively. The benchmark data are shown as open circles. Note that the 
standard- and revised- algorithm curves are indistinguishable in (b) at h — 2, 4 and 
6 fs. 
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